######################################################################
### Fig 9: Lower Tariffs on Differentiated Products
######################################################################
rm(list=ls())

## setting the working directory to the current folder
setwd(getwd())

## libraries
library(arm)
library(Matrix)
library(lme4)
library(mvtnorm)

data <- read.csv("./hs8_sigma_ntb_BG.csv",
                 stringsAsFactors=F)

## id variable for jt
a <- paste(as.character(data$naics), as.character(data$year), sep="_")
data$id_jt <- a
data$id_jt <- as.factor(data$id_jt)
data$year <- as.factor(data$year)
data$naics <- as.factor(data$naics)

rmna <- which(is.na(data$sigma))
sub <- data[-rmna,]

## subset only manufacturing industry
nrow(sub)
sub <- sub[sub$manuf==1,]
nrow(sub)

D <- subset(sub, select=c("AdValorem", "year", "naics","sigma", "id_jt","low.sigma", "high.sigma", "med.sigma",
                   "naics.emp", "naics.vadd", "naics.tfp4", "naics.pay",
                   "naics.energy", "ad", "cvd", "hs8.imp.value", "hs8.n.imp.cty"))

## remove NAs
D <- na.omit(D)

lme <- lmer(AdValorem ~ low.sigma+high.sigma + log(naics.emp)+ log(naics.vadd)+
            log(naics.tfp4)+ log(naics.pay)+ log(naics.energy)+
            ad + cvd + hs8.imp.value + hs8.n.imp.cty +
            (1+low.sigma+high.sigma|year) + (1|id_jt) + (1|naics),
            data=D)

summary(lme)


#########################################################################
### quasi-bayesian simulation for moving from high sigma to low sigma

### fixed effects
fe.mean <- fixef(lme)
fe.var <- as.matrix(vcov(lme))

### random effects
## 5381 naics-year
re_idjt.mean <- ranef(lme)$id_jt[,1]
re_idjt.var <- as.vector(attr(ranef(lme,postVar=T)[[1]], "postVar"))
## 373 naics
re_naics.mean <- ranef(lme)$naics[,1]
re_naics.var <- as.vector(attr(ranef(lme,postVar=T)[[2]], "postVar"))
## 15 years
re_time.mean <- ranef(lme)$year
re_time.var <- attr(ranef(lme,postVar=T)[[3]], "postVar")

nsim <- 1000
re_idjt.draw <- rep(NA, length(re_idjt.mean))
names(re_idjt.draw) <- rownames(ranef(lme)$id_jt)

re_naics.draw <- rep(NA, length(re_naics.mean))
names(re_naics.draw) <- rownames(ranef(lme)$naics)

re_time.draw <- matrix(NA, nrow=nrow(re_time.mean), ncol=ncol(re_time.mean))
rownames(re_time.draw) <- rownames(re_time.mean)

effect <- matrix(NA, ncol=nsim, nrow=nrow(re_time.mean))
rownames(effect) <- rownames(re_time.draw)

D$year <- as.character(D$year)

#########################################################################
### simulation
for(i in 1:nsim){
  if (i %% 10==0){
    print(i)
    flush.console()
  }
  fe.draw <- rmvnorm(1, fe.mean, fe.var, method="svd")

  for(j in 1:length(re_idjt.mean)) {
    re_idjt.draw[j] <- rnorm(1, re_idjt.mean[j], sqrt(re_idjt.var[j]))
  }

  for(j in 1:length(re_naics.mean)) {
    re_naics.draw[j] <- rnorm(1, re_naics.mean[j], sqrt(re_naics.var[j]))
  }

  
  for(j in 1:nrow(re_time.mean)) {
    re_time.draw[j,] <- rmvnorm(1, mean=as.numeric(re_time.mean[j,]),
                               sigma=re_time.var[,,j], method="svd")
  }
  
  for(y in unique(as.character(D$year))){

    X.1 <- X.0 <- model.matrix(lme)[D$year==y,]

    X.1[,2] <- rep(1, nrow(X.1)) # treatment: low sigma
    X.1[,3] <- rep(0, nrow(X.1)) # treatment: low sigma

    X.0[,2] <- rep(0, nrow(X.1)) # treatment: low sigma
    X.0[,3] <- rep(1, nrow(X.1)) # treatment: low sigma

    eta.1 <- as.numeric(fe.draw %*% t(X.1)) +
      as.numeric(re_idjt.draw[match(D$id_jt[D$year==y],names(re_idjt.draw))]) +
        as.numeric(re_naics.draw[match(D$naics[D$year==y],names(re_naics.draw))]) + 
          as.numeric(re_time.draw[rownames(re_time.draw)==y,] %*% c(1,1,0))
    
    eta.0 <- as.numeric(fe.draw %*% t(X.0)) +
      as.numeric(re_idjt.draw[match(D$id_jt[D$year==y],names(re_idjt.draw))]) +
        as.numeric(re_naics.draw[match(D$naics[D$year==y],names(re_naics.draw))]) + 
          as.numeric(re_time.draw[rownames(re_time.draw)==y,] %*% c(1,0,1))

    
    effect[rownames(effect)==y,i] <- mean(eta.1) - mean(eta.0)
  }
}


## Plot estimated coefficient
effect.year <- apply(effect, 1, mean)
effect.year.95up <- apply(effect, 1, quantile, .975)
effect.year.95lo <- apply(effect, 1, quantile, .025)

every <- c(effect.year.95up, effect.year.95lo)
years <- as.vector(rownames(ranef(lme)$year))

pdf(file="./Figure9.pdf",
    width=10, height=8)
par(oma=c(0,5,0,1))
plot(years, effect.year, pch=16,
     xlab="", ylab="", xaxt='n', yaxt='n',
     ylim=c(-0.5,0.2))

for(i in 1:length(years)){
  lines(c(years[i],years[i]), c(effect.year.95up[i],effect.year.95lo[i]))
}

abline(h=0, col="red")
abline(v=1995, col="gray", lty=2)

text(1992, -0.4, "Uruguay Round", cex=1.5)
arrows(1990, -0.35, 1995, -0.35 , length = 0.25, angle = 30,
       code = 2, cex=.9, col="gray")


years <- sort(unique(c(years, "1994")))
years <- seq(1990,2005,3)
axis(1, at=years, labels=rep("", length(years)))
for (i in 1:length(years)){
  ## xlabel
  xloc.i <- as.numeric(years[i])+.5
  text(xloc.i, min(effect.year.95lo)-.06, adj=1, labels=years[i], xpd=T, cex=1.7)
}

mtext(side=2, "Effect of Product Differentiation on Ad-valorem Tariff\n(Percentage Point Change)",
      line=4.9, cex=1.7)
tmp <- c((min(effect.year.95lo)), (max(effect.year.95up)))
tmp <- round(tmp, digits=2)
tmp[2] <- 0.2
tmp2 <- round(seq(tmp[1], tmp[2], 0.1), digits=2)
axis(2, at=tmp2, label=rep("",length(tmp2)))
for (i in 1:length(tmp2)){
  ## xlabel
  yloc.i <- tmp2[i]
  text(1988.9, yloc.i, adj=1, labels=tmp2[i], xpd=T, cex=1.7)
}


dev.off()


